/[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.1 by mocchiut, Fri May 19 13:15:56 2006 UTC revision 1.3 by pam-fi, Thu Oct 26 16:22:37 2006 UTC
# Line 12  Line 12 
12  ************************************************************************    ************************************************************************  
13    
14    
15        SUBROUTINE MINI_2(ISTEP,IFAIL)        SUBROUTINE MINI2(ISTEP,IFAIL,IPRINT)
16                
17        IMPLICIT DOUBLE PRECISION (A-H,O-Z)        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
18    
19        include 'commontracker.f' !tracker general common        include 'commontracker.f' !tracker general common
20        include 'common_mini_2.f' !common for the tracking procedure        include 'common_mini_2.f' !common for the tracking procedure
21    
22        logical DEBUG  c      logical DEBUG
23        common/dbg/DEBUG  c      common/dbg/DEBUG
24    
25        parameter (inf=1.e8)      !just a huge number...        parameter (inf=1.e8)      !just a huge number...
26  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 31  c     the plane ordering is reversed in Line 31  c     the plane ordering is reversed in
31  c     ordering, but they maintain their Z coordinates. so plane number 1 is  c     ordering, but they maintain their Z coordinates. so plane number 1 is
32  c     the first one that a particle meets, and its Z coordinate is > 0  c     the first one that a particle meets, and its Z coordinate is > 0
33  c------------------------------------------------------------------------  c------------------------------------------------------------------------
34        DATA ZINI/23.5/           !z coordinate of the reference plane        DATA ZINI/23.5/  !!! ***PP*** to be changed !z coordinate of the reference plane
35    
36        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/100/        !maximum number of steps in the chi^2 minimization        DATA ISTEPMAX/100/        !maximum number of steps in the chi^2 minimization
# Line 48  c      DATA ALMIN/-inf,-inf,-inf,-inf,-0 Line 48  c      DATA ALMIN/-inf,-inf,-inf,-inf,-0
48    
49    
50        DIMENSION DAL(5)                    !increment of vector alfa        DIMENSION DAL(5)                    !increment of vector alfa
51          DIMENSION CHI2DD_R(4,4),CHI2D_R(4) !hessiano e gradiente di chi2
52    
53        INTEGER IFLAG                    INTEGER IFLAG            
54  c--------------------------------------------------------  c--------------------------------------------------------
# Line 59  c Line 60  c
60  c     NB: the two metods gives equivalent results BUT  c     NB: the two metods gives equivalent results BUT
61  c     method 2 is faster!!  c     method 2 is faster!!
62  c--------------------------------------------------------  c--------------------------------------------------------
63        DATA IFLAG/2/                    DATA IFLAG/2/  
64            
65          LOGICAL TRKDEBUG
66          COMMON/TRKD/TRKDEBUG
67    
68          IF(IPRINT.EQ.1) THEN
69             TRKDEBUG = .TRUE.
70          ELSE
71             TRKDEBUG = .FALSE.
72          ENDIF        
73    
74  *     ----------------------------------------------------------  *     ----------------------------------------------------------
75  *     define ALTOL(5) ---> tolerances on state vector  *     define ALTOL(5) ---> tolerances on state vector
# Line 79  c     deflection error (see PDG) Line 89  c     deflection error (see PDG)
89  *      *    
90        ISTEP=0                   !num. steps to minimize chi^2        ISTEP=0                   !num. steps to minimize chi^2
91        JFAIL=0                   !error flag        JFAIL=0                   !error flag
       CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives  
       IF(JFAIL.NE.0) THEN  
          IFAIL=1  
          if(DEBUG)  
      $        PRINT *,'mini_2 ===> error on CHISQ computation !!!'  
          RETURN  
       ENDIF  
92  *      *    
93  *     -----------------------  *     -----------------------
94  *     START MINIMIZATION LOOP  *     START MINIMIZATION LOOP
95  *     -----------------------  *     -----------------------
96   10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!   10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!
97        CHI2_P=CHI2  
98          CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives
99          IF(JFAIL.NE.0) THEN
100             IFAIL=1
101             CHI2=-9999.
102             if(TRKDEBUG)
103         $        PRINT *,'*** ERROR in mini *** wrong CHISQ'
104             RETURN
105          ENDIF
106          
107    *      CHI2_P=CHI2
108                
109  c      print*,'@@@@@ ',istep,' - ',al  c      print*,'@@@@@ ',istep,' - ',al
110    
111        cost=1e-7        COST=1e-7
112        DO I=1,5        DO I=1,5
113           DO J=1,5           DO J=1,5
114              CHI2DD(I,J)=CHI2DD(I,J)*COST              CHI2DD(I,J)=CHI2DD(I,J)*COST
115           ENDDO           ENDDO
116           CHI2D(I)=CHI2D(I)*COST           CHI2D(I)=CHI2D(I)*COST
117        ENDDO              ENDDO      
118    
119          IF(PFIXED.EQ.0.) THEN
120    
121  *------------------------------------------------------------*  *------------------------------------------------------------*
122  *     track fitting with FREE deflection  *     track fitting with FREE deflection
123  *------------------------------------------------------------*  *------------------------------------------------------------*
124        CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
125        IF(IFA.NE.0) THEN         !not positive-defined                 IF(IFA.NE.0) THEN      !not positive-defined      
126           if(DEBUG)then              if(TRKDEBUG)then
127              PRINT *,                 PRINT *,
128       $     'MINI_HOUGH ==> '//       $              '*** ERROR in mini ***'//
129       $     '** ERROR ** on matrix inversion (not positive-defined)!!!'       $              'on matrix inversion (not pos-def)'
130       $     ,DET       $              ,DET
131           endif              endif
132           IFAIL=1              IF(CHI2.EQ.0) CHI2=-9999.
133           RETURN                          IF(CHI2.GT.0) CHI2=-CHI2
134        ENDIF              IFAIL=1
135        CALL DSFINV(5,CHI2DD,5)   !CHI2DD matrix inversion                  RETURN            
136             ENDIF
137             CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
138  *     *******************************************  *     *******************************************
139  *     find new value of AL-pha                 !*  *     find new value of AL-pha                
140  *                                              !*  *                                              
141        DO I=1,5                                 !*           DO I=1,5                              
142           DAL(I)=0.                             !*              DAL(I)=0.          
143           DO J=1,5                              !*              DO J=1,5            
144              DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J) !*                 DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)
145              COV(I,J)=CHI2DD(I,J)               !*                 COV(I,J)=2.*COST*CHI2DD(I,J)
146           ENDDO                                 !*              ENDDO              
147        ENDDO                                    !*           ENDDO                  
148        DO I=1,5                                 !*           DO I=1,5              
149           AL(I)=AL(I)+DAL(I)                    !*              AL(I)=AL(I)+DAL(I)  
150        ENDDO                                    !*           ENDDO                  
151    *------------------------------------------------------------*
152    *     track fitting with FIXED deflection
153    *------------------------------------------------------------*
154          ELSE
155             AL(5)=1./PFIXED
156             DO I=1,4
157                CHI2D_R(I)=CHI2D(I)
158                DO J=1,4
159                   CHI2DD_R(I,J)=CHI2DD(I,J)
160                ENDDO
161             ENDDO
162             CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
163             IF(IFA.NE.0) THEN
164                if(TRKDEBUG)then
165                   PRINT *,
166         $              '*** ERROR in mini ***'//
167         $              'on matrix inversion (not pos-def)'
168         $              ,DET
169                endif
170                IF(CHI2.EQ.0) CHI2=-9999.
171                IF(CHI2.GT.0) CHI2=-CHI2
172                IFAIL=1
173                RETURN            
174             ENDIF
175             CALL DSFINV(4,CHI2DD_R,4)
176             DO I=1,4
177                DAL(I)=0.
178                DO J=1,4
179                   DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J)
180                   COV(I,J)=2.*COST*CHI2DD_R(I,J)
181                ENDDO
182             ENDDO
183             DAL(5)=0.
184             DO I=1,4
185                AL(I)=AL(I)+DAL(I)
186             ENDDO
187          ENDIF
188    *------------------------------------------------------------*
189    *     ----------------------------------------------------   *
190    *------------------------------------------------------------*
191  *     *******************************************  *     *******************************************
192  *     check parameter bounds:  *     check parameter bounds:
193        DO I=1,5        DO I=1,5
194           IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN           IF(AL(I).GT.ALMAX(I).OR.AL(I).LT.ALMIN(I))THEN
195              if(DEBUG)then              if(TRKDEBUG)then
196                 PRINT*,' **WARNING**  '                 PRINT*,' *** WARNING in mini ***  '
197                 PRINT*,'MINI_2 ==> AL(',I,') out of range'                 PRINT*,'MINI_2 ==> AL(',I,') out of range'
198                 PRINT*,'         value: ',AL(I),                 PRINT*,'         value: ',AL(I),
199       $              '  limits: ',ALMIN(I),ALMAX(I)                     $              '  limits: ',ALMIN(I),ALMAX(I)              
200                 print*,'istep ',istep                 print*,'istep ',istep
201              endif              endif
202                IF(CHI2.EQ.0) CHI2=-9999.
203                IF(CHI2.GT.0) CHI2=-CHI2
204              IFAIL=1              IFAIL=1
205              RETURN              RETURN
206           ENDIF           ENDIF
207        ENDDO              ENDDO      
208  *     new estimate of chi^2:  
       JFAIL=0                   !error flag  
       CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives  
       IF(JFAIL.NE.0) THEN  
          IFAIL=1  
          if(DEBUG)  
      $        PRINT *,'mini_2:  ===> error on CHISQ computation !!!'  
          RETURN  
       ENDIF  
209  *     check number of steps:  *     check number of steps:
210        IF(ISTEP.gt.ISTEPMAX) then        IF(ISTEP.gt.ISTEPMAX) then
211           IFAIL=1           IFAIL=1
212           if(DEBUG)           if(TRKDEBUG)
213       $        PRINT *,'mini_2: WARNING ===> ISTEP.GT.ISTEPMAX=',ISTEPMAX       $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
214         $        ISTEPMAX
215           goto 11           goto 11
216        endif        endif
217  *     ---------------------------------------------  *     ---------------------------------------------
# Line 171  ',istep,' - ',al Line 224  ',istep,' - ',al
224           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
225        ENDDO        ENDDO
226    
227    *     new estimate of chi^2:
228          JFAIL=0                   !error flag
229          CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives
230          IF(JFAIL.NE.0) THEN
231             IFAIL=1
232             if(TRKDEBUG)THEN
233                CHI2=-9999.
234             if(TRKDEBUG)
235         $        PRINT *,'*** ERROR in mini *** wrong CHISQ'
236             ENDIF
237             RETURN
238          ENDIF
239          COST=1e-7
240          DO I=1,5
241             DO J=1,5
242                CHI2DD(I,J)=CHI2DD(I,J)*COST
243             ENDDO
244             CHI2D(I)=CHI2D(I)*COST
245          ENDDO      
246          IF(PFIXED.EQ.0.) THEN
247             CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
248             IF(IFA.NE.0) THEN      !not positive-defined      
249                if(TRKDEBUG)then
250                   PRINT *,
251         $              '*** ERROR in mini ***'//
252         $              'on matrix inversion (not pos-def)'
253         $              ,DET
254                endif
255                IF(CHI2.EQ.0) CHI2=-9999.
256                IF(CHI2.GT.0) CHI2=-CHI2
257                IFAIL=1
258                RETURN            
259             ENDIF
260             CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
261             DO I=1,5                              
262                DAL(I)=0.          
263                DO J=1,5            
264                   COV(I,J)=2.*COST*CHI2DD(I,J)
265                ENDDO              
266             ENDDO                  
267          ELSE
268             DO I=1,4
269                CHI2D_R(I)=CHI2D(I)
270                DO J=1,4
271                   CHI2DD_R(I,J)=CHI2DD(I,J)
272                ENDDO
273             ENDDO
274             CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
275             IF(IFA.NE.0) THEN
276                if(TRKDEBUG)then
277                   PRINT *,
278         $              '*** ERROR in mini ***'//
279         $              'on matrix inversion (not pos-def)'
280         $              ,DET
281                endif
282                IF(CHI2.EQ.0) CHI2=-9999.
283                IF(CHI2.GT.0) CHI2=-CHI2
284                IFAIL=1
285                RETURN            
286             ENDIF
287             CALL DSFINV(4,CHI2DD_R,4)
288             DO I=1,4
289                DAL(I)=0.
290                DO J=1,4
291                   COV(I,J)=2.*COST*CHI2DD_R(I,J)
292                ENDDO
293             ENDDO
294          ENDIF
295    *****************************
296    
297  *     ------------------------------------  *     ------------------------------------
298  *     Number of Degree Of Freedom  *     Number of Degree Of Freedom
# Line 180  ',istep,' - ',al Line 302  ',istep,' - ',al
302       $        +int(xgood(ip))       $        +int(xgood(ip))
303       $        +int(ygood(ip))       $        +int(ygood(ip))
304        enddo        enddo
305        ndof=ndof-5                                  if(pfixed.eq.0.) ndof=ndof-5 ! ***PP***                          
306          if(pfixed.ne.0.) ndof=ndof-4 ! ***PP***
307          if(ndof.le.0.) then
308             ndof = 1
309             if(TRKDEBUG)
310         $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'
311          endif
312  *     ------------------------------------  *     ------------------------------------
313  *     Reduced chi^2  *     Reduced chi^2
314        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
315    
         
316   11   CONTINUE         11   CONTINUE      
317    
318   101  CONTINUE        NSTEP=ISTEP ! ***PP***
   
 c      print*,'END MINI'  
319    
320        RETURN        RETURN
321        END                END        
# Line 217  c      print*,'END MINI' Line 342  c      print*,'END MINI'
342        DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes)        DIMENSION XV2(nplanes),YV2(nplanes),XV1(nplanes),YV1(nplanes)
343       $     ,XV0(nplanes),YV0(nplanes)       $     ,XV0(nplanes),YV0(nplanes)
344        DIMENSION AL_P(5)        DIMENSION AL_P(5)
345    
346          LOGICAL TRKDEBUG
347          COMMON/TRKD/TRKDEBUG
348  *      *    
349  *     chi^2 computation  *     chi^2 computation
350  *      *    
# Line 226  c      print*,'END MINI' Line 354  c      print*,'END MINI'
354        JFAIL=0                   !error flag        JFAIL=0                   !error flag
355        CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes        CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes
356        IF(JFAIL.NE.0) THEN        IF(JFAIL.NE.0) THEN
357           PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!'           IF(TRKDEBUG)
358         $        PRINT *,'CHISQ ==> error from trk routine POSXYZ !!'
359           IFAIL=1           IFAIL=1
360           RETURN           RETURN
361        ENDIF        ENDIF
# Line 297  c$$$      ENDDO Line 426  c$$$      ENDDO
426              JFAIL=0              JFAIL=0
427              CALL POSXYZ(AL_P,JFAIL)              CALL POSXYZ(AL_P,JFAIL)
428              IF(JFAIL.NE.0) THEN              IF(JFAIL.NE.0) THEN
429                 PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!'                 IF(TRKDEBUG)
430    *23456789012345678901234567890123456789012345678901234567890123456789012
431         $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
432                 IFAIL=1                 IFAIL=1
433                 RETURN                 RETURN
434              ENDIF              ENDIF
# Line 309  c$$$      ENDDO Line 440  c$$$      ENDDO
440              JFAIL=0              JFAIL=0
441              CALL POSXYZ(AL_P,JFAIL)              CALL POSXYZ(AL_P,JFAIL)
442              IF(JFAIL.NE.0) THEN              IF(JFAIL.NE.0) THEN
443                 PRINT *,'CHISQ ==> error from tracking routine POSXYZ !!'                 IF(TRKDEBUG)
444         $              PRINT *,'CHISQ ==> error from trk routine POSXYZ'
445                 IFAIL=1                 IFAIL=1
446                 RETURN                 RETURN
447              ENDIF              ENDIF
# Line 340  c$$$      ENDDO Line 472  c$$$      ENDDO
472                            
473              COSTHE=DSQRT(1.-AL(3)**2)              COSTHE=DSQRT(1.-AL(3)**2)
474              IF(COSTHE.EQ.0.) THEN              IF(COSTHE.EQ.0.) THEN
475                 PRINT *,'=== WARNING ===> COSTHE=0'                 IF(TRKDEBUG)PRINT *,'=== WARNING ===> COSTHE=0'
476                 STOP                 IFAIL=1
477                   RETURN
478              ENDIF              ENDIF
479                            
480              DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3              DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3
# Line 425  c$$$      ENDDO Line 558  c$$$      ENDDO
558                
559        include 'commontracker.f' !tracker general common        include 'commontracker.f' !tracker general common
560        include 'common_mini_2.f' !common for the tracking procedure        include 'common_mini_2.f' !common for the tracking procedure
561    
562          LOGICAL TRKDEBUG
563          COMMON/TRKD/TRKDEBUG
564  c        c      
565        DIMENSION AL_P(5)        DIMENSION AL_P(5)
566  *      *    
# Line 454  c       Line 590  c      
590           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)           CALL GRKUTA(CHARGE,STEP,VECT,VOUT)
591           IF(VOUT(3).GT.VECT(3)) THEN           IF(VOUT(3).GT.VECT(3)) THEN
592              IFAIL=1              IFAIL=1
593              PRINT *,'posxy (grkuta): WARNING ===> backward track!!'              if(TRKDEBUG)
594              print*,'charge',charge       $      PRINT *,'posxy (grkuta): WARNING ===> backward track!!'
595              print*,'vect',vect              if(.TRUE.)print*,'charge',charge
596              print*,'vout',vout              if(.TRUE.)print*,'vect',vect
597              print*,'step',step              if(.TRUE.)print*,'vout',vout
598                if(.TRUE.)print*,'step',step
599              RETURN              RETURN
600           ENDIF           ENDIF
601           Z=VOUT(3)           Z=VOUT(3)
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23