/[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.4 by pam-fi, Tue Nov 7 15:55:11 2006 UTC revision 1.18 by pam-fi, Fri Jun 1 15:01:19 2007 UTC
# Line 36  c--------------------------------------- Line 36  c---------------------------------------
36  c      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/120/        !maximum number of steps in the chi^2 minimization        DATA ISTEPMAX/100/        !maximum number of steps in the chi^2 minimization
40        DATA TOLL/1.d-8/          !tolerance in reaching the next plane during        DATA TOLL/1.d-8/          !tolerance in reaching the next plane during
41  *                               !the tracking procedure  *                               !the tracking procedure
42        DATA STEPMAX/100./        !maximum number of steps in the trackin gprocess        DATA STEPMAX/100./        !maximum number of steps in the trackin gprocess
43    
44    c      DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
45    c      DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
46        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components        DATA ALMAX/dinf,dinf,1.,dinf,dinf/ !limits on alpha vector components
47        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"        DATA ALMIN/-dinf,-dinf,-1.,-dinf,-dinf/ !"
48    
49        DIMENSION DAL(5)                    !increment of vector alfa  c$$$      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
51    
52  c     elena--------  c     elena--------
# Line 65  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        LOGICAL TRKDEBUG,TRKVERBOSE,STUDENT,FIRSTSTEPS,FIRSTSTUDENT
71        COMMON/TRKD/TRKDEBUG,TRKVERBOSE        COMMON/TRKD/TRKDEBUG,TRKVERBOSE
72    
73          DIMENSION AL0(5)
74          LOGICAL SUCCESS_NEW,SUCCESS_OLD
75    *
76    *     define kind of minimization (0x=chi2+gaussian or 1x=likelihood+student)
77    *
78          STUDENT = .false.
79          FIRSTSTEPS = .true.
80          FIRSTSTUDENT = .true.
81          IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true.
82    
83        IF(IPRINT.EQ.1) THEN        IF(IPRINT.EQ.1) THEN
84           TRKVERBOSE = .TRUE.           TRKVERBOSE = .TRUE.
85           TRKDEBUG   = .FALSE.           TRKDEBUG   = .FALSE.
# Line 103  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE Line 115  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
115  *     ----------------------------------------------------------  *     ----------------------------------------------------------
116  *     changed in order to evaluate energy-dependent  *     changed in order to evaluate energy-dependent
117  *     tolerances on all 5 parameters  *     tolerances on all 5 parameters
118        FACT=100.                  !scale factor to define tolerance on alfa  cPP      FACT=1.0e10                  !scale factor to define tolerance on alfa
119  c     deflection error (see PDG)  c     deflection error (see PDG)
120        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.))
121        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 119  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451 Line 131  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
131  *      *    
132        ISTEP=0                   !num. steps to minimize chi^2        ISTEP=0                   !num. steps to minimize chi^2
133        JFAIL=0                   !error flag        JFAIL=0                   !error flag
134          CHI2=0
135    
136  c     elena--------        if(TRKDEBUG) print*,'guess: ',al
       CHI2   = 0  
 c     elena--------  
137        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
138    
139  *      *    
# Line 131  c     elena-------- Line 142  c     elena--------
142  *     -----------------------  *     -----------------------
143   10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!   10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!
144    
145        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives  * -------------------------------
146        IF(JFAIL.NE.0) THEN  * **** Chi2+gaussian minimization
147           IFAIL=1  * -------------------------------
148           CHI2=-9999.  
149           if(TRKVERBOSE)        IF(.NOT.STUDENT.OR.FIRSTSTEPS) THEN
150       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'  
151           RETURN           IF(ISTEP.GE.3) FIRSTSTEPS = .false.
152        ENDIF  
153                   CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
154        COST=1e-7           IF(JFAIL.NE.0) THEN
155        DO I=1,5              IFAIL=1
156           DO J=1,5              CHI2=-9999.
157              CHI2DD(I,J)=CHI2DD(I,J)*COST              if(TRKVERBOSE)
158         $           PRINT *,'*** ERROR in mini *** wrong CHISQ'
159                RETURN
160             ENDIF
161            
162    c      COST=1e-5
163             COST=1.
164             DO I=1,5
165                IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2
166             ENDDO
167             DO I=1,5
168                DO J=1,5
169                   CHI2DD(I,J)=CHI2DD(I,J)*COST
170                ENDDO
171    c$$$            CHI2D(I)=CHI2D(I)*COST
172           ENDDO           ENDDO
          CHI2D(I)=CHI2D(I)*COST  
       ENDDO        
173    
174        IF(PFIXED.EQ.0.) THEN           IF(PFIXED.EQ.0.) THEN
175    
176  *------------------------------------------------------------*  *------------------------------------------------------------*
177  *     track fitting with FREE deflection  *     track fitting with FREE deflection
178  *------------------------------------------------------------*  *------------------------------------------------------------*
179           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant              CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
180           IF(IFA.NE.0) THEN      !not positive-defined                    IF(IFA.NE.0) THEN   !not positive-defined      
181              if(TRKVERBOSE)then                 if(TRKVERBOSE)then
182                 PRINT *,                    PRINT *,
183       $              '*** ERROR in mini ***'//       $                 '*** ERROR in mini ***'//
184       $              'on matrix inversion (not pos-def)'       $                 'on matrix inversion (not pos-def)'
185       $              ,DET       $                 ,DET
186              endif                 endif
187              IF(CHI2.EQ.0) CHI2=-9999.                 IF(CHI2.EQ.0) CHI2=-9999.
188              IF(CHI2.GT.0) CHI2=-CHI2                 IF(CHI2.GT.0) CHI2=-CHI2
189              IFAIL=1                 IFAIL=1
190              RETURN                             RETURN            
191           ENDIF              ENDIF
192           CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion                  CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
193  *     *******************************************  *     *******************************************
194  *     find new value of AL-pha                  *     find new value of AL-pha                
195  *     *******************************************                                          *     *******************************************                                        
196           DO I=1,5                                            DO I=1,5                              
197              DAL(I)=0.                           DAL(I)=0.          
198              DO J=1,5                             DO J=1,5            
199                 DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)                    DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)  *COST
200                 COV(I,J)=2.*COST*CHI2DD(I,J)                    COV(I,J)=2.*COST*CHI2DD(I,J)
201              ENDDO                               ENDDO  
202           ENDDO                                ENDDO                  
203           DO I=1,5                            DO I=1,5              
204              AL(I)=AL(I)+DAL(I)                   AL(I)=AL(I)+DAL(I)  
205           ENDDO                                ENDDO                  
206  *------------------------------------------------------------*  *------------------------------------------------------------*
207  *     track fitting with FIXED deflection  *     track fitting with FIXED deflection
208  *------------------------------------------------------------*  *------------------------------------------------------------*
209        ELSE           ELSE
210           AL(5)=1./PFIXED              AL(5)=1./PFIXED
211           DO I=1,4              DO I=1,4
212              CHI2D_R(I)=CHI2D(I)                 CHI2D_R(I)=CHI2D(I)
213              DO J=1,4                 DO J=1,4
214                 CHI2DD_R(I,J)=CHI2DD(I,J)                    CHI2DD_R(I,J)=CHI2DD(I,J)
215                   ENDDO
216              ENDDO              ENDDO
217           ENDDO              CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
218           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)              IF(IFA.NE.0) THEN
219           IF(IFA.NE.0) THEN                 if(TRKVERBOSE)then
220              if(TRKVERBOSE)then                    PRINT *,
221                 PRINT *,       $                 '*** ERROR in mini ***'//
222       $              '*** ERROR in mini ***'//       $                 'on matrix inversion (not pos-def)'
223       $              'on matrix inversion (not pos-def)'       $                 ,DET
224       $              ,DET                 endif
225              endif                 IF(CHI2.EQ.0) CHI2=-9999.
226              IF(CHI2.EQ.0) CHI2=-9999.                 IF(CHI2.GT.0) CHI2=-CHI2
227              IF(CHI2.GT.0) CHI2=-CHI2                 IFAIL=1
228              IFAIL=1                 RETURN            
229              RETURN                          ENDIF
230           ENDIF              CALL DSFINV(4,CHI2DD_R,4)
          CALL DSFINV(4,CHI2DD_R,4)  
231  *     *******************************************  *     *******************************************
232  *     find new value of AL-pha                  *     find new value of AL-pha                
233  *     *******************************************                                          *     *******************************************                                        
234           DO I=1,4              DO I=1,4
235              DAL(I)=0.                 DAL(I)=0.
236              DO J=1,4                 DO J=1,4
237                 DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J)                    DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J)  *COST
238                 COV(I,J)=2.*COST*CHI2DD_R(I,J)                    COV(I,J)=2.*COST*CHI2DD_R(I,J)
239                   ENDDO
240                ENDDO
241                DAL(5)=0.
242                DO I=1,4
243                   AL(I)=AL(I)+DAL(I)
244                ENDDO
245             ENDIF
246          
247             if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
248    
249    c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
250    c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
251    
252    
253          ENDIF
254    
255    * -------------------------------
256    * **** Likelihood+Student minimization
257    * -------------------------------
258    
259          IF(STUDENT.AND.(.NOT.FIRSTSTEPS)) THEN
260    
261             IF(FIRSTSTUDENT) THEN
262                FIRSTSTUDENT = .false.
263                ISTEP = 1
264             ENDIF
265    
266             CALL CHISQSTT(1,JFAIL)
267             DO I=1,5                              
268                DAL(I)=0.          
269                DO J=1,5            
270                   DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)
271                ENDDO  
272             ENDDO                  
273            
274             DO I=1,5
275                DO j=1,5
276                   COV(I,J) = 2.*CHI2DD(I,J)
277              ENDDO              ENDDO
278           ENDDO           ENDDO
279           DAL(5)=0.          
280           DO I=1,4           CHI2TOLL = 1.E-3
281              AL(I)=AL(I)+DAL(I)           ALPHA = 3.0
282             BETA = -0.4
283             E=1.
284             EA=1.
285             EB=1.
286             EC=1.
287             FA=1.
288             FB=1.
289             FC=1.
290             SUCCESS_OLD = .FALSE.
291             SUCCESS_NEW = .FALSE.
292    
293             CALL CHISQSTT(0,JFAIL)
294    c$$$         PRINT*,CHI2
295             CHI2_NEW = CHI2
296             FC = CHI2
297             EC = 0.
298    
299     100     CONTINUE
300             DO I=1,5
301                AL0(I)=AL(I)
302             ENDDO
303             DO I=1,5
304                AL(I)=AL(I)+E*DAL(I)
305           ENDDO           ENDDO
306             CALL CHISQSTT(0,JFAIL)
307             CHI2_OLD = CHI2_NEW
308             CHI2_NEW = CHI2
309             FA = FB
310             FB = FC
311             FC = CHI2
312             EA = EB
313             EB = EC
314             EC = E
315    
316    c$$$         PRINT*,E,CHI2_NEW
317    
318             IF(CHI2_NEW.LE.CHI2_OLD) THEN ! success
319                IF(DABS(CHI2_NEW-CHI2_OLD).LT.CHI2TOLL) GOTO 101
320                SUCCESS_OLD = SUCCESS_NEW
321                SUCCESS_NEW = .TRUE.
322                E = E*ALPHA
323             ELSE                   ! failure
324                SUCCESS_OLD = SUCCESS_NEW
325                SUCCESS_NEW = .FALSE.
326                CHI2_NEW = CHI2_OLD
327                DO I=1,5
328                   AL(I)=AL0(I)
329                ENDDO
330                IF(SUCCESS_OLD) THEN
331                   DENOM = (EB-EA)*(FB-FC) - (EB-EC)*(FB-FA)
332                   IF(DENOM.NE.0.) THEN
333                      E = EB - 0.5*( (EB-EA)**2*(FB-FC)
334         $                 - (EB-EC)**2*(FB-FA) ) / DENOM
335                   ELSE
336                      E = BETA*E
337                   ENDIF
338                ELSE
339                   E = BETA*E
340                ENDIF
341    c$$$            E = BETA*E
342             ENDIF
343             GOTO 100
344    
345     101     CONTINUE
346    
347             DO I=1,5
348                DAL(I)=E*DAL(I)
349             ENDDO
350    
351    c$$$         print*,' '
352    c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
353    c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
354    c$$$         print*,'==== CHI2 ===='
355    c$$$         print*,chi2
356    c$$$         print*,'==== CHI2d ===='
357    c$$$         print*,(chi2d(i),i=1,5)
358    c$$$         print*,'==== CHI2dd ===='
359    c$$$         do j=1,5
360    c$$$            print*,(chi2dd(j,i),i=1,5)
361    c$$$         enddo
362    c$$$         print*,'================'
363    c$$$         print*,' '
364            
365    *========= FIN QUI =============
366    
367        ENDIF        ENDIF
368    
369        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)          
370    
371    
372    
373  *------------------------------------------------------------*  *------------------------------------------------------------*
374  *     ----------------------------------------------------   *  *     ----------------------------------------------------   *
# Line 247  c     elena-------- Line 394  c     elena--------
394  *     check number of steps:  *     check number of steps:
395  *------------------------------------------------------------*  *------------------------------------------------------------*
396        IF(ISTEP.ge.ISTEPMAX) then        IF(ISTEP.ge.ISTEPMAX) then
397           IFAIL=1  c$$$         IFAIL=1
398           if(TRKVERBOSE)  c$$$         if(TRKVERBOSE)
399       $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',  c$$$     $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
400       $        ISTEPMAX  c$$$     $        ISTEPMAX
401           goto 11           goto 11
402        endif        endif
403  *------------------------------------------------------------*  *------------------------------------------------------------*
# Line 266  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DE Line 413  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DE
413        ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51        ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51
414        ALTOL(4) = ALTOL(3)                ALTOL(4) = ALTOL(3)        
415                
416    c$$$      print*,' -- ',(DAL(I),ALTOL(I),' - ',i=1,5) !>>>> new step!
417    
418  *---- check tolerances:  *---- check tolerances:
419  c$$$      DO I=1,5  c$$$      DO I=1,5
420  c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!  c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!
421  c$$$      ENDDO  c$$$      ENDDO
422  c$$$      print*,'chi2 -- ',DCHI2  c$$$      print*,'chi2 -- ',DCHI2
423    
424          IF(ISTEP.LT.ISTEPMIN) GOTO 10 ! ***PP***
425        DO I=1,5        DO I=1,5
426           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
427        ENDDO        ENDDO
428    
429  *     new estimate of chi^2:  *****************************
430        JFAIL=0                   !error flag  *     final estimate of chi^2
431        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives  *****************************
432        IF(JFAIL.NE.0) THEN  
433           IFAIL=1  * -------------------------------
434           if(TRKVERBOSE)THEN  * **** Chi2+gaussian minimization
435              CHI2=-9999.  * -------------------------------
436           if(TRKVERBOSE)  
437       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'        IF(.NOT.STUDENT) THEN
438           ENDIF  
439           RETURN           JFAIL=0                !error flag
440        ENDIF           CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
441        COST=1e-7           IF(JFAIL.NE.0) THEN
       DO I=1,5  
          DO J=1,5  
             CHI2DD(I,J)=CHI2DD(I,J)*COST  
          ENDDO  
          CHI2D(I)=CHI2D(I)*COST  
       ENDDO        
       IF(PFIXED.EQ.0.) THEN  
          CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant  
          IF(IFA.NE.0) THEN      !not positive-defined        
             if(TRKVERBOSE)then  
                PRINT *,  
      $              '*** ERROR in mini ***'//  
      $              'on matrix inversion (not pos-def)'  
      $              ,DET  
             endif  
             IF(CHI2.EQ.0) CHI2=-9999.  
             IF(CHI2.GT.0) CHI2=-CHI2  
442              IFAIL=1              IFAIL=1
443              RETURN                          if(TRKVERBOSE)THEN
444                   CHI2=-9999.
445                   if(TRKVERBOSE)
446         $              PRINT *,'*** ERROR in mini *** wrong CHISQ'
447                ENDIF
448                RETURN
449           ENDIF           ENDIF
450           CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion      c      COST=1e-7
451           DO I=1,5                                         COST=1.
452              DAL(I)=0.                     DO I=1,5
453              DO J=1,5                          IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2
454                 COV(I,J)=2.*COST*CHI2DD(I,J)           ENDDO
455              ENDDO                         DO I=1,5
456           ENDDO                                DO J=1,5
457        ELSE                 CHI2DD(I,J)=CHI2DD(I,J)*COST
          DO I=1,4  
             CHI2D_R(I)=CHI2D(I)  
             DO J=1,4  
                CHI2DD_R(I,J)=CHI2DD(I,J)  
458              ENDDO              ENDDO
459           ENDDO           ENDDO
460           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)           IF(PFIXED.EQ.0.) THEN
461           IF(IFA.NE.0) THEN              CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
462              if(TRKVERBOSE)then              IF(IFA.NE.0) THEN   !not positive-defined      
463                 PRINT *,                 if(TRKVERBOSE)then
464       $              '*** ERROR in mini ***'//                    PRINT *,
465       $              'on matrix inversion (not pos-def)'       $                 '*** ERROR in mini ***'//
466       $              ,DET       $                 'on matrix inversion (not pos-def)'
467              endif       $                 ,DET
468              IF(CHI2.EQ.0) CHI2=-9999.                 endif
469              IF(CHI2.GT.0) CHI2=-CHI2                 IF(CHI2.EQ.0) CHI2=-9999.
470              IFAIL=1                 IF(CHI2.GT.0) CHI2=-CHI2
471              RETURN                             IFAIL=1
472                   RETURN            
473                ENDIF
474                CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
475                DO I=1,5                              
476    c$$$               DAL(I)=0.          
477                   DO J=1,5            
478                      COV(I,J)=2.*COST*CHI2DD(I,J)
479                   ENDDO              
480                ENDDO                  
481             ELSE
482                DO I=1,4
483                   CHI2D_R(I)=CHI2D(I)
484                   DO J=1,4
485                      CHI2DD_R(I,J)=CHI2DD(I,J)
486                   ENDDO
487                ENDDO
488                CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
489                IF(IFA.NE.0) THEN
490                   if(TRKVERBOSE)then
491                      PRINT *,
492         $                 '*** ERROR in mini ***'//
493         $                 'on matrix inversion (not pos-def)'
494         $                 ,DET
495                   endif
496                   IF(CHI2.EQ.0) CHI2=-9999.
497                   IF(CHI2.GT.0) CHI2=-CHI2
498                   IFAIL=1
499                   RETURN            
500                ENDIF
501                CALL DSFINV(4,CHI2DD_R,4)
502                DO I=1,4
503    c$$$               DAL(I)=0.
504                   DO J=1,4
505                      COV(I,J)=2.*COST*CHI2DD_R(I,J)
506                   ENDDO
507                ENDDO
508           ENDIF           ENDIF
509           CALL DSFINV(4,CHI2DD_R,4)  
510           DO I=1,4        ENDIF
511              DAL(I)=0.  
512              DO J=1,4  * -------------------------------
513                 COV(I,J)=2.*COST*CHI2DD_R(I,J)  * **** Likelihood+student minimization
514    * -------------------------------
515    
516          IF(STUDENT) THEN
517             CALL CHISQSTT(1,JFAIL)
518             DO I=1,5
519                DO j=1,5
520                   COV(I,J) = 2.*CHI2DD(I,J)
521              ENDDO              ENDDO
522           ENDDO           ENDDO
523        ENDIF        ENDIF
524    
525  *****************************  *****************************
526    
527  *     ------------------------------------  *     ------------------------------------
# Line 362  c$$$      print*,'chi2 -- ',DCHI2 Line 540  c$$$      print*,'chi2 -- ',DCHI2
540       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'
541        endif        endif
542    
       if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)  
   
543  *     ------------------------------------  *     ------------------------------------
544  *     Reduced chi^2  *     Reduced chi^2
545        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
# Line 372  c      print*,'mini2: chi2 ',chi2 Line 548  c      print*,'mini2: chi2 ',chi2
548    
549   11   CONTINUE         11   CONTINUE      
550    
551          if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)
552    
553        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
554    
555    c$$$      print*,'>>>>> NSTEP = ',NSTEP
556    
557        RETURN        RETURN
558        END                END        
559    
# Line 439  c$$$      ENDDO Line 619  c$$$      ENDDO
619  *     measured position of the cluster.  *     measured position of the cluster.
620  *     ---------------------------------------------------------  *     ---------------------------------------------------------
621        CHI2=0.        CHI2=0.
622        DO I=1,nplanes        DO I=1,nplanes        
623           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
624              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))
625              ALFA = XM_A(I) - BETA * YM_A(I)              ALFA = XM_A(I) - BETA * YM_A(I)
# Line 466  c$$$      ENDDO Line 646  c$$$      ENDDO
646       +                                       *( XGOOD(I)*(1-YGOOD(I)) )       +                                       *( XGOOD(I)*(1-YGOOD(I)) )
647       +        +((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
648       +                                       *( (1-XGOOD(I))*YGOOD(I) )       +                                       *( (1-XGOOD(I))*YGOOD(I) )
649    c$$$         print*,(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
650    c$$$         print*,(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
651    c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
652    c$$$     +        *( XGOOD(I)*(1-YGOOD(I)) )
653    c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
654    c$$$     +        *( (1-XGOOD(I))*YGOOD(I) )
655    c$$$         print*,XV(I),XM(I),XGOOD(I)
656    c$$$         print*,YV(I),YM(I),YGOOD(I)
657        ENDDO        ENDDO
658    c$$$      print*,'CHISQ ',chi2
659  *     ------------------------------------------------  *     ------------------------------------------------
660  *      *    
661  *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)      *     calculation of derivatives (dX/dAL_fa and dY/dAL_fa)    
# Line 586  c$$$      ENDDO Line 775  c$$$      ENDDO
775        RETURN        RETURN
776        END        END
777    
778    ******************************************************************************    
779    *    
780    *     routine to compute Likelihodd+Student and its derivatives
781    *    
782    *     (modified in respect to the previous one in order to include
783    *     single clusters. In this case the residual is evaluated by
784    *     calculating the distance between the track intersection and the
785    *     segment AB associated to the single cluster)
786    *
787    ******************************************************************************
788    
789          SUBROUTINE CHISQSTT(IFLAG,JFAIL)
790    
791          IMPLICIT DOUBLE PRECISION (A-H,O-Z)
792    
793          include 'commontracker.f' !tracker general common
794          include 'common_mini_2.f' !common for the tracking procedure
795    
796          LOGICAL TRKDEBUG,TRKVERBOSE
797          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
798    
799          DIMENSION AL_P(5)
800          DIMENSION VECTEMP(5)
801    c$$$      DIMENSION U(5) ! BFGS
802    
803          DO I=1,5
804             AL_P(I)=AL(I)
805          ENDDO
806          JFAIL=0                   !error flag
807          CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes
808          IF(JFAIL.NE.0) THEN
809             IF(TRKVERBOSE)
810         $        PRINT *,'CHISQSTT ==> error from trk routine POSXYZ !!'
811             IFAIL=1
812             RETURN
813          ENDIF
814    
815          DO I=1,nplanes
816             DXDAL(I,1)=1.
817             DYDAL(I,1)=0.
818             DXDAL(I,2)=0.
819             DYDAL(I,2)=1.
820             COSTHE=DSQRT(1.-AL(3)**2)
821             IF(COSTHE.EQ.0.) THEN
822                IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
823                IFAIL=1
824                RETURN
825             ENDIF
826             DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3
827             DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3
828             DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE
829             DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE
830             IF(AL(5).NE.0.) THEN      
831                DXDAL(I,5)=
832         +           (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I))
833         +           *DCOS(AL(4))))/AL(5)
834                DYDAL(I,5)=
835         +           (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I))
836         +           *DSIN(AL(4))))/AL(5)
837             ELSE
838                DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 )
839                DYDAL(I,5)=0.
840             ENDIF
841          ENDDO
842    
843          IF(IFLAG.EQ.0) THEN ! function calulation
844             CHI2=0.
845             DO I=1,nplanes
846                IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
847                   BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
848                   ALFA = XM_A(I) - BETA * YM_A(I)
849                   YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
850                   if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
851         $              YM(I)=dmin1(YM_A(I),YM_B(I))
852                   if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
853         $              YM(I)=dmax1(YM_A(I),YM_B(I))
854                   XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
855                ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
856                   BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
857                   ALFA = YM_A(I) - BETA * XM_A(I)
858                   XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
859                   if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
860         $              XM(I)=dmin1(XM_A(I),XM_B(I))
861                   if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
862         $              XM(I)=dmax1(XM_A(I),XM_B(I))
863                   YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
864                ENDIF
865                TERMX = DLOG( (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)/
866         $           (TAILX(I)*RESX(I)**2) )
867                TERMY = DLOG( (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)/
868         $           (TAILY(I)*RESY(I)**2) )
869                CHI2=CHI2
870         $           +(TAILX(I)+1.0)*TERMX   *( XGOOD(I) )
871         $           +(TAILY(I)+1.0)*TERMY   *( YGOOD(I) )
872             ENDDO
873          ENDIF
874    
875          IF(IFLAG.EQ.1) THEN ! derivative calulation
876             DO I=1,5
877                CHI2DOLD(I)=CHI2D(I)
878             ENDDO
879             DO J=1,5
880                CHI2D(J)=0.      
881                DO I=1,nplanes
882                   CHI2D(J)=CHI2D(J)
883         $              +2.*(TAILX(I)+1.0)*(XV(I)-XM(I))/
884         $              (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)*
885         $              DXDAL(I,J)   *XGOOD(I)
886         $              +2.*(TAILY(I)+1.0)*(YV(I)-YM(I))/
887         $              (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)*
888         $              DYDAL(I,J)   *YGOOD(I)
889                ENDDO
890             ENDDO
891             DO K=1,5
892                VECTEMP(K)=0.
893                DO M=1,5
894                   VECTEMP(K) = VECTEMP(K) +
895         $              COV(K,M)/2.*(CHI2D(M)-CHI2DOLD(M))
896                ENDDO
897             ENDDO
898             DOWN1 = 0.
899             DO K=1,5
900                DOWN1 = DOWN1 + DAL(K)*(CHI2D(K)-CHI2DOLD(K))
901             ENDDO
902             IF(DOWN1.EQ.0.) THEN
903                PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN1 = 0'
904                IFAIL=1
905                RETURN
906             ENDIF
907             DOWN2 = 0.
908             DO K=1,5
909                DO M=1,5
910                   DOWN2 = DOWN2 + (CHI2D(K)-CHI2DOLD(K))*VECTEMP(K)
911                ENDDO
912             ENDDO
913             IF(DOWN2.EQ.0.) THEN
914                PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN2 = 0'
915                IFAIL=1
916                RETURN
917             ENDIF
918    c$$$         DO K=1,5 ! BFGS
919    c$$$            U(K) = DAL(K)/DOWN1 - VECTEMP(K)/DOWN2
920    c$$$         ENDDO
921             DO I=1,5
922                DO J=1,5
923                   CHI2DD(I,J) = COV(I,J)/2.
924         $              +DAL(I)*DAL(J)/DOWN1
925         $              -VECTEMP(I)*VECTEMP(J)/DOWN2
926    c$$$     $              +DOWN2*U(I)*U(J) ! BFGS
927                ENDDO
928             ENDDO
929          ENDIF
930    
931          RETURN
932          END
933          
934  *****************************************************************  *****************************************************************
935  *  *
936  *     Routine to compute the track intersection points  *     Routine to compute the track intersection points
# Line 625  c      COMMON/TRKD/TRKVERBOSE Line 969  c      COMMON/TRKD/TRKVERBOSE
969  c        c      
970        DIMENSION AL_P(5)        DIMENSION AL_P(5)
971  *      *    
972        DO I=1,nplanes  cpp      DO I=1,nplanes
973           ZV(I)=ZM(I)            !  cpp         ZV(I)=ZM(I)            !
974        ENDDO        cpp      ENDDO      
975  *      *    
976  *     set parameters for GRKUTA  *     set parameters for GRKUTA
977  *      *    
# Line 641  c       Line 985  c      
985        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
986        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))
987        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
988    
989    c$$$      print*,'POSXY (prima) ',vout
990    
991        DO I=1,nplanes        DO I=1,nplanes
992           step=vout(3)-zv(i)  cpp         step=vout(3)-zv(i)
993             step=vout(3)-zm(i)
994   10      DO J=1,7   10      DO J=1,7
995              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
996              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
# Line 674  c$$$            if(.TRUE.)print*,'step', Line 1022  c$$$            if(.TRUE.)print*,'step',
1022              GOTO 11              GOTO 11
1023           ENDIF           ENDIF
1024    
1025    
1026  *     -----------------------------------------------  *     -----------------------------------------------
1027  *        evaluate track coordinates  *        evaluate track coordinates
1028   100     XV(I)=VOUT(1)   100     XV(I)=VOUT(1)
# Line 683  c$$$            if(.TRUE.)print*,'step', Line 1032  c$$$            if(.TRUE.)print*,'step',
1032           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
1033  *     -----------------------------------------------  *     -----------------------------------------------
1034    
1035             IF(TRACKMODE.EQ.1) THEN
1036    *     -----------------------------------------------
1037    *        change of energy by bremsstrahlung for electrons
1038                VOUT(7) = VOUT(7) * 0.997 !0.9968
1039    *     -----------------------------------------------
1040             ENDIF
1041    
1042        ENDDO        ENDDO
1043    
1044    c$$$      print*,'POSXY (dopo) ',vout
1045    
1046    
1047        RETURN        RETURN
1048        END        END
1049    
# Line 725  c         ZM_B(IP) = 0 Line 1084  c         ZM_B(IP) = 0
1084           RESY(IP) = 1000.       !12.d-4           RESY(IP) = 1000.       !12.d-4
1085           XGOOD(IP) = 0           XGOOD(IP) = 0
1086           YGOOD(IP) = 0           YGOOD(IP) = 0
1087             DEDXTRK_X(IP) = 0
1088             DEDXTRK_Y(IP) = 0
1089             AXV(IP) = 0
1090             AYV(IP) = 0
1091             XV(IP) = -100
1092             YV(IP) = -100
1093        enddo        enddo
1094    
1095        return        return
# Line 795  c      IMPLICIT DOUBLE PRECISION (A-H,O- Line 1160  c      IMPLICIT DOUBLE PRECISION (A-H,O-
1160              ZP(NP)=ZM(I)              ZP(NP)=ZM(I)
1161           ENDIF           ENDIF
1162        ENDDO        ENDDO
1163          IFLAG=0                   !no debug mode
1164        CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)        CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)
1165  c      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI  
1166    c$$$      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG
1167    c$$$      print*,' XP ',(xp(i),i=1,np)
1168    c$$$      print*,' ZP ',(zp(i),i=1,np)
1169    c$$$      print*,' AP ',(ap(i),i=1,np)
1170    c$$$      print*,' XP ',(rp(i),i=1,np)
1171    
1172        IF(IFLAG.NE.0)GOTO 10 !straigth fit        IF(IFLAG.NE.0)GOTO 10 !straigth fit
1173    c      if(CHI.gt.100)GOTO 10 !straigth fit
1174        ARG = RADIUS**2-(ZINI-ZC)**2        ARG = RADIUS**2-(ZINI-ZC)**2
1175        IF(ARG.LT.0)GOTO 10       !straigth fit        IF(ARG.LT.0)GOTO 10       !straigth fit
1176        DC = SQRT(ARG)              DC = SQRT(ARG)      
# Line 806  c      print*,' circle: ',XC,ZC,RADIUS,' Line 1179  c      print*,' circle: ',XC,ZC,RADIUS,'
1179        AX = -(ZINI-ZC)/DC        AX = -(ZINI-ZC)/DC
1180        DEF=100./(RADIUS*0.3*0.43)        DEF=100./(RADIUS*0.3*0.43)
1181        IF(XC.GT.0)DEF=-DEF        IF(XC.GT.0)DEF=-DEF
1182          
1183    
1184    
1185          IF(ABS(X0).GT.30)THEN
1186    c$$$         PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS
1187    c$$$     $     ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF
1188             GOTO 10       !straigth fit
1189          ENDIF
1190        GOTO 20                   !guess is ok        GOTO 20                   !guess is ok
1191    
1192  *     ----------------------------------------  *     ----------------------------------------
# Line 847  c      print*,' circle: ',XC,ZC,RADIUS,' Line 1228  c      print*,' circle: ',XC,ZC,RADIUS,'
1228        AL(2) = Y0        AL(2) = Y0
1229        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
1230        AL(3) = tath/sqrt(1+tath**2)        AL(3) = tath/sqrt(1+tath**2)
1231        IF(AX.NE.0)THEN  c$$$      IF(AX.NE.0)THEN
1232           AL(4)= atan(AY/AX)  c$$$         AL(4)= atan(AY/AX)
1233        ELSE  c$$$      ELSE
1234           AL(4) = acos(-1.)/2  c$$$         AL(4) = acos(-1.)/2
1235           IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)  c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
1236    c$$$  ENDIF
1237    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
1238    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
1239    
1240    c$$$      AL(4) =  0.
1241    c$$$      IF(AX.NE.0.AND.AY.NE.0)THEN
1242    c$$$         AL(4)= atan(AY/AX)
1243    c$$$      ELSEIF(AY.EQ.0)THEN
1244    c$$$         AL(4) = 0.
1245    c$$$         IF(AX.LT.0)AL(4) = AL(4)+acos(-1.)        
1246    c$$$      ELSEIF(AX.EQ.0)THEN
1247    c$$$         AL(4) = acos(-1.)/2
1248    c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
1249    c$$$      ENDIF
1250    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
1251    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
1252    
1253    c$$$      AL(4)=0.
1254    c$$$      IF( AX.NE.0.OR.AY.NE.0. ) THEN
1255    c$$$         AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
1256    c$$$         IF(AX.LT.0.) AL(4) = ACOS(-1.0)-AL(4)
1257    c$$$      ENDIF
1258    
1259          AL(4)=0.
1260          IF( AX.NE.0.OR.AY.NE.0. ) THEN
1261             AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
1262             IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4)
1263             IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4)
1264        ENDIF        ENDIF
1265        IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)        IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0)
1266        AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs        IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0)
1267    
1268        AL(5) = DEF        AL(5) = DEF
1269    
1270  c      print*,' guess: ',(al(i),i=1,5)  c      print*,' guess: ',(al(i),i=1,5)
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.23