/[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.17 by pam-fi, Thu May 24 13:29:09 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
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          IF(MOD(INT(TRACKMODE/10),10).EQ.1) STUDENT = .true.
80    
81        IF(IPRINT.EQ.1) THEN        IF(IPRINT.EQ.1) THEN
82           TRKVERBOSE = .TRUE.           TRKVERBOSE = .TRUE.
83           TRKDEBUG   = .FALSE.           TRKDEBUG   = .FALSE.
# Line 103  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE Line 113  c      COMMON/TRKD/TRKDEBUG,TRKVERBOSE
113  *     ----------------------------------------------------------  *     ----------------------------------------------------------
114  *     changed in order to evaluate energy-dependent  *     changed in order to evaluate energy-dependent
115  *     tolerances on all 5 parameters  *     tolerances on all 5 parameters
116        FACT=100.                  !scale factor to define tolerance on alfa  cPP      FACT=1.0e10                  !scale factor to define tolerance on alfa
117  c     deflection error (see PDG)  c     deflection error (see PDG)
118        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.))
119        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 129  c$$$      DELETA2 = 0.016/0.3/0.4/0.4451
129  *      *    
130        ISTEP=0                   !num. steps to minimize chi^2        ISTEP=0                   !num. steps to minimize chi^2
131        JFAIL=0                   !error flag        JFAIL=0                   !error flag
132          CHI2=0
133    
134  c     elena--------        if(TRKDEBUG) print*,'guess: ',al
       CHI2   = 0  
 c     elena--------  
135        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
136    
137  *      *    
# Line 131  c     elena-------- Line 140  c     elena--------
140  *     -----------------------  *     -----------------------
141   10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!   10   ISTEP=ISTEP+1             !<<<<<<<<<<<<<< NEW STEP !!
142    
143        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives  * -------------------------------
144        IF(JFAIL.NE.0) THEN  * **** Chi2+gaussian minimization
145           IFAIL=1  * -------------------------------
146           CHI2=-9999.  
147           if(TRKVERBOSE)        IF(.NOT.STUDENT) THEN
148       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'  
149           RETURN           CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
150        ENDIF           IF(JFAIL.NE.0) THEN
151                      IFAIL=1
152        COST=1e-7              CHI2=-9999.
153        DO I=1,5              if(TRKVERBOSE)
154           DO J=1,5       $           PRINT *,'*** ERROR in mini *** wrong CHISQ'
155              CHI2DD(I,J)=CHI2DD(I,J)*COST              RETURN
156             ENDIF
157            
158    c      COST=1e-5
159             COST=1.
160             DO I=1,5
161                IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2
162             ENDDO
163             DO I=1,5
164                DO J=1,5
165                   CHI2DD(I,J)=CHI2DD(I,J)*COST
166                ENDDO
167    c$$$            CHI2D(I)=CHI2D(I)*COST
168           ENDDO           ENDDO
          CHI2D(I)=CHI2D(I)*COST  
       ENDDO        
169    
170        IF(PFIXED.EQ.0.) THEN           IF(PFIXED.EQ.0.) THEN
171    
172  *------------------------------------------------------------*  *------------------------------------------------------------*
173  *     track fitting with FREE deflection  *     track fitting with FREE deflection
174  *------------------------------------------------------------*  *------------------------------------------------------------*
175           CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant              CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
176           IF(IFA.NE.0) THEN      !not positive-defined                    IF(IFA.NE.0) THEN   !not positive-defined      
177              if(TRKVERBOSE)then                 if(TRKVERBOSE)then
178                 PRINT *,                    PRINT *,
179       $              '*** ERROR in mini ***'//       $                 '*** ERROR in mini ***'//
180       $              'on matrix inversion (not pos-def)'       $                 'on matrix inversion (not pos-def)'
181       $              ,DET       $                 ,DET
182              endif                 endif
183              IF(CHI2.EQ.0) CHI2=-9999.                 IF(CHI2.EQ.0) CHI2=-9999.
184              IF(CHI2.GT.0) CHI2=-CHI2                 IF(CHI2.GT.0) CHI2=-CHI2
185              IFAIL=1                 IFAIL=1
186              RETURN                             RETURN            
187           ENDIF              ENDIF
188           CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion                  CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
189  *     *******************************************  *     *******************************************
190  *     find new value of AL-pha                  *     find new value of AL-pha                
191  *     *******************************************                                          *     *******************************************                                        
192           DO I=1,5                                            DO I=1,5                              
193              DAL(I)=0.                           DAL(I)=0.          
194              DO J=1,5                             DO J=1,5            
195                 DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)                    DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)  *COST
196                 COV(I,J)=2.*COST*CHI2DD(I,J)                    COV(I,J)=2.*COST*CHI2DD(I,J)
197              ENDDO                               ENDDO  
198           ENDDO                                ENDDO                  
199           DO I=1,5                            DO I=1,5              
200              AL(I)=AL(I)+DAL(I)                   AL(I)=AL(I)+DAL(I)  
201           ENDDO                                ENDDO                  
202  *------------------------------------------------------------*  *------------------------------------------------------------*
203  *     track fitting with FIXED deflection  *     track fitting with FIXED deflection
204  *------------------------------------------------------------*  *------------------------------------------------------------*
205        ELSE           ELSE
206           AL(5)=1./PFIXED              AL(5)=1./PFIXED
207           DO I=1,4              DO I=1,4
208              CHI2D_R(I)=CHI2D(I)                 CHI2D_R(I)=CHI2D(I)
209              DO J=1,4                 DO J=1,4
210                 CHI2DD_R(I,J)=CHI2DD(I,J)                    CHI2DD_R(I,J)=CHI2DD(I,J)
211                   ENDDO
212              ENDDO              ENDDO
213           ENDDO              CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
214           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)              IF(IFA.NE.0) THEN
215           IF(IFA.NE.0) THEN                 if(TRKVERBOSE)then
216              if(TRKVERBOSE)then                    PRINT *,
217                 PRINT *,       $                 '*** ERROR in mini ***'//
218       $              '*** ERROR in mini ***'//       $                 'on matrix inversion (not pos-def)'
219       $              'on matrix inversion (not pos-def)'       $                 ,DET
220       $              ,DET                 endif
221              endif                 IF(CHI2.EQ.0) CHI2=-9999.
222              IF(CHI2.EQ.0) CHI2=-9999.                 IF(CHI2.GT.0) CHI2=-CHI2
223              IF(CHI2.GT.0) CHI2=-CHI2                 IFAIL=1
224              IFAIL=1                 RETURN            
225              RETURN                          ENDIF
226           ENDIF              CALL DSFINV(4,CHI2DD_R,4)
          CALL DSFINV(4,CHI2DD_R,4)  
227  *     *******************************************  *     *******************************************
228  *     find new value of AL-pha                  *     find new value of AL-pha                
229  *     *******************************************                                          *     *******************************************                                        
230           DO I=1,4              DO I=1,4
231              DAL(I)=0.                 DAL(I)=0.
232              DO J=1,4                 DO J=1,4
233                 DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J)                    DAL(I)=DAL(I)-CHI2DD_R(I,J)*CHI2D_R(J)  *COST
234                 COV(I,J)=2.*COST*CHI2DD_R(I,J)                    COV(I,J)=2.*COST*CHI2DD_R(I,J)
235                   ENDDO
236                ENDDO
237                DAL(5)=0.
238                DO I=1,4
239                   AL(I)=AL(I)+DAL(I)
240              ENDDO              ENDDO
241             ENDIF
242          
243             if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)
244    
245    c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
246    c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
247    
248    
249          ENDIF
250    
251    * -------------------------------
252    * **** Likelihood+Student minimization
253    * -------------------------------
254    
255          IF(STUDENT) THEN
256             CALL CHISQSTT(1,JFAIL)
257             DO I=1,5                              
258                DAL(I)=0.          
259                DO J=1,5            
260                   DAL(I)=DAL(I)-CHI2DD(I,J)*CHI2D(J)
261                ENDDO  
262             ENDDO                  
263            
264             DO I=1,5
265                DO j=1,5
266                   COV(I,J) = 2.*CHI2DD(I,J)
267                ENDDO
268             ENDDO
269            
270             CHI2TOLL = 1.E-3
271             ALPHA = 3.0
272             BETA = -0.4
273             E=1.
274             EA=1.
275             EB=1.
276             EC=1.
277             FA=1.
278             FB=1.
279             FC=1.
280             SUCCESS_OLD = .FALSE.
281             SUCCESS_NEW = .FALSE.
282    
283             CALL CHISQSTT(0,JFAIL)
284    c$$$         PRINT*,CHI2
285             CHI2_NEW = CHI2
286             FC = CHI2
287             EC = 0.
288    
289     100     CONTINUE
290             DO I=1,5
291                AL0(I)=AL(I)
292           ENDDO           ENDDO
293           DAL(5)=0.           DO I=1,5
294           DO I=1,4              AL(I)=AL(I)+E*DAL(I)
             AL(I)=AL(I)+DAL(I)  
295           ENDDO           ENDDO
296             CALL CHISQSTT(0,JFAIL)
297             CHI2_OLD = CHI2_NEW
298             CHI2_NEW = CHI2
299             FA = FB
300             FB = FC
301             FC = CHI2
302             EA = EB
303             EB = EC
304             EC = E
305    
306    c$$$         PRINT*,E,CHI2_NEW
307    
308             IF(CHI2_NEW.LE.CHI2_OLD) THEN ! success
309                IF(DABS(CHI2_NEW-CHI2_OLD).LT.CHI2TOLL) GOTO 101
310                SUCCESS_OLD = SUCCESS_NEW
311                SUCCESS_NEW = .TRUE.
312                E = E*ALPHA
313             ELSE                   ! failure
314                SUCCESS_OLD = SUCCESS_NEW
315                SUCCESS_NEW = .FALSE.
316                CHI2_NEW = CHI2_OLD
317                DO I=1,5
318                   AL(I)=AL0(I)
319                ENDDO
320                IF(SUCCESS_OLD) THEN
321                   DENOM = (EB-EA)*(FB-FC) - (EB-EC)*(FB-FA)
322                   IF(DENOM.NE.0.) THEN
323                      E = EB - 0.5*( (EB-EA)**2*(FB-FC)
324         $                 - (EB-EC)**2*(FB-FA) ) / DENOM
325                   ELSE
326                      E = BETA*E
327                   ENDIF
328                ELSE
329                   E = BETA*E
330                ENDIF
331    c$$$            E = BETA*E
332             ENDIF
333             GOTO 100
334    
335     101     CONTINUE
336    
337             DO I=1,5
338                DAL(I)=E*DAL(I)
339             ENDDO
340    
341    c$$$         print*,' '
342    c$$$         PRINT*,'DAL ',(DAL(K),K=1,5)
343    c$$$         PRINT*,'CHI2DOLD ',(CHI2DOLD(K),K=1,5)
344    c$$$         print*,'==== CHI2 ===='
345    c$$$         print*,chi2
346    c$$$         print*,'==== CHI2d ===='
347    c$$$         print*,(chi2d(i),i=1,5)
348    c$$$         print*,'==== CHI2dd ===='
349    c$$$         do j=1,5
350    c$$$            print*,(chi2dd(j,i),i=1,5)
351    c$$$         enddo
352    c$$$         print*,'================'
353    c$$$         print*,' '
354            
355    *========= FIN QUI =============
356    
357        ENDIF        ENDIF
358    
359        if(TRKDEBUG) print*,'mini2: step ',istep,chi2,1./AL(5)          
360    
361    
362    
363  *------------------------------------------------------------*  *------------------------------------------------------------*
364  *     ----------------------------------------------------   *  *     ----------------------------------------------------   *
# Line 247  c     elena-------- Line 384  c     elena--------
384  *     check number of steps:  *     check number of steps:
385  *------------------------------------------------------------*  *------------------------------------------------------------*
386        IF(ISTEP.ge.ISTEPMAX) then        IF(ISTEP.ge.ISTEPMAX) then
387           IFAIL=1  c$$$         IFAIL=1
388           if(TRKVERBOSE)  c$$$         if(TRKVERBOSE)
389       $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',  c$$$     $        PRINT *,'*** WARNING in mini *** ISTEP.GT.ISTEPMAX=',
390       $        ISTEPMAX  c$$$     $        ISTEPMAX
391           goto 11           goto 11
392        endif        endif
393  *------------------------------------------------------------*  *------------------------------------------------------------*
# Line 266  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DE Line 403  c$$$      ALTOL(5) = DSQRT(DELETA1**2+DE
403        ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51        ALTOL(3) = DSQRT(ALTOL(1)**2+ALTOL(2)**2)/44.51
404        ALTOL(4) = ALTOL(3)                ALTOL(4) = ALTOL(3)        
405                
406    c$$$      print*,' -- ',(DAL(I),ALTOL(I),' - ',i=1,5) !>>>> new step!
407    
408  *---- check tolerances:  *---- check tolerances:
409  c$$$      DO I=1,5  c$$$      DO I=1,5
410  c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!  c$$$         if(TRKVERBOSE)print*,i,' -- ',DAL(I),ALTOL(I) !>>>> new step!
411  c$$$      ENDDO  c$$$      ENDDO
412  c$$$      print*,'chi2 -- ',DCHI2  c$$$      print*,'chi2 -- ',DCHI2
413    
414          IF(ISTEP.LT.ISTEPMIN) GOTO 10 ! ***PP***
415        DO I=1,5        DO I=1,5
416           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!           IF(ABS(DAL(I)).GT.ALTOL(I))GOTO 10 !>>>> new step!
417        ENDDO        ENDDO
418    
419  *     new estimate of chi^2:  *****************************
420        JFAIL=0                   !error flag  *     final estimate of chi^2
421        CALL CHISQ(IFLAG,JFAIL)   !chi^2 and its derivatives  *****************************
422        IF(JFAIL.NE.0) THEN  
423           IFAIL=1  * -------------------------------
424           if(TRKVERBOSE)THEN  * **** Chi2+gaussian minimization
425              CHI2=-9999.  * -------------------------------
426           if(TRKVERBOSE)  
427       $        PRINT *,'*** ERROR in mini *** wrong CHISQ'        IF(.NOT.STUDENT) THEN
428           ENDIF  
429           RETURN           JFAIL=0                !error flag
430        ENDIF           CALL CHISQ(IFLAG,JFAIL) !chi^2 and its derivatives
431        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  
432              IFAIL=1              IFAIL=1
433              RETURN                          if(TRKVERBOSE)THEN
434                   CHI2=-9999.
435                   if(TRKVERBOSE)
436         $              PRINT *,'*** ERROR in mini *** wrong CHISQ'
437                ENDIF
438                RETURN
439           ENDIF           ENDIF
440           CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion      c      COST=1e-7
441           DO I=1,5                                         COST=1.
442              DAL(I)=0.                     DO I=1,5
443              DO J=1,5                          IF(CHI2DD(I,I).NE.0.)COST=COST/DABS(CHI2DD(I,I))**0.2
444                 COV(I,J)=2.*COST*CHI2DD(I,J)           ENDDO
445              ENDDO                         DO I=1,5
446           ENDDO                                DO J=1,5
447        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)  
448              ENDDO              ENDDO
449           ENDDO           ENDDO
450           CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)           IF(PFIXED.EQ.0.) THEN
451           IF(IFA.NE.0) THEN              CALL DSFACT(5,CHI2DD,5,IFA,DET,JFA) !CHI2DD matrix determinant
452              if(TRKVERBOSE)then              IF(IFA.NE.0) THEN   !not positive-defined      
453                 PRINT *,                 if(TRKVERBOSE)then
454       $              '*** ERROR in mini ***'//                    PRINT *,
455       $              'on matrix inversion (not pos-def)'       $                 '*** ERROR in mini ***'//
456       $              ,DET       $                 'on matrix inversion (not pos-def)'
457              endif       $                 ,DET
458              IF(CHI2.EQ.0) CHI2=-9999.                 endif
459              IF(CHI2.GT.0) CHI2=-CHI2                 IF(CHI2.EQ.0) CHI2=-9999.
460              IFAIL=1                 IF(CHI2.GT.0) CHI2=-CHI2
461              RETURN                             IFAIL=1
462                   RETURN            
463                ENDIF
464                CALL DSFINV(5,CHI2DD,5) !CHI2DD matrix inversion    
465                DO I=1,5                              
466    c$$$               DAL(I)=0.          
467                   DO J=1,5            
468                      COV(I,J)=2.*COST*CHI2DD(I,J)
469                   ENDDO              
470                ENDDO                  
471             ELSE
472                DO I=1,4
473                   CHI2D_R(I)=CHI2D(I)
474                   DO J=1,4
475                      CHI2DD_R(I,J)=CHI2DD(I,J)
476                   ENDDO
477                ENDDO
478                CALL DSFACT(4,CHI2DD_R,4,IFA,DET,JFA)
479                IF(IFA.NE.0) THEN
480                   if(TRKVERBOSE)then
481                      PRINT *,
482         $                 '*** ERROR in mini ***'//
483         $                 'on matrix inversion (not pos-def)'
484         $                 ,DET
485                   endif
486                   IF(CHI2.EQ.0) CHI2=-9999.
487                   IF(CHI2.GT.0) CHI2=-CHI2
488                   IFAIL=1
489                   RETURN            
490                ENDIF
491                CALL DSFINV(4,CHI2DD_R,4)
492                DO I=1,4
493    c$$$               DAL(I)=0.
494                   DO J=1,4
495                      COV(I,J)=2.*COST*CHI2DD_R(I,J)
496                   ENDDO
497                ENDDO
498           ENDIF           ENDIF
499           CALL DSFINV(4,CHI2DD_R,4)  
500           DO I=1,4        ENDIF
501              DAL(I)=0.  
502              DO J=1,4  * -------------------------------
503                 COV(I,J)=2.*COST*CHI2DD_R(I,J)  * **** Likelihood+student minimization
504    * -------------------------------
505    
506          IF(STUDENT) THEN
507             CALL CHISQSTT(1,JFAIL)
508             DO I=1,5
509                DO j=1,5
510                   COV(I,J) = 2.*CHI2DD(I,J)
511              ENDDO              ENDDO
512           ENDDO           ENDDO
513        ENDIF        ENDIF
514    
515  *****************************  *****************************
516    
517  *     ------------------------------------  *     ------------------------------------
# Line 362  c$$$      print*,'chi2 -- ',DCHI2 Line 530  c$$$      print*,'chi2 -- ',DCHI2
530       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'       $        print*,'*** WARNING *** in mini n.dof = 0 (set to 1)'
531        endif        endif
532    
       if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)  
   
533  *     ------------------------------------  *     ------------------------------------
534  *     Reduced chi^2  *     Reduced chi^2
535        CHI2 = CHI2/dble(ndof)        CHI2 = CHI2/dble(ndof)
# Line 372  c      print*,'mini2: chi2 ',chi2 Line 538  c      print*,'mini2: chi2 ',chi2
538    
539   11   CONTINUE         11   CONTINUE      
540    
541          if(TRKDEBUG) print*,'mini2: -ok- ',istep,chi2,1./AL(5)
542    
543        NSTEP=ISTEP ! ***PP***        NSTEP=ISTEP ! ***PP***
544    
545    c$$$      print*,'>>>>> NSTEP = ',NSTEP
546    
547        RETURN        RETURN
548        END                END        
549    
# Line 439  c$$$      ENDDO Line 609  c$$$      ENDDO
609  *     measured position of the cluster.  *     measured position of the cluster.
610  *     ---------------------------------------------------------  *     ---------------------------------------------------------
611        CHI2=0.        CHI2=0.
612        DO I=1,nplanes        DO I=1,nplanes        
613           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
614              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))
615              ALFA = XM_A(I) - BETA * YM_A(I)              ALFA = XM_A(I) - BETA * YM_A(I)
# Line 466  c$$$      ENDDO Line 636  c$$$      ENDDO
636       +                                       *( XGOOD(I)*(1-YGOOD(I)) )       +                                       *( XGOOD(I)*(1-YGOOD(I)) )
637       +        +((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
638       +                                       *( (1-XGOOD(I))*YGOOD(I) )       +                                       *( (1-XGOOD(I))*YGOOD(I) )
639    c$$$         print*,(XV(I)-XM(I))**2/RESX(i)**2   *( XGOOD(I)*YGOOD(I) )
640    c$$$         print*,(YV(I)-YM(I))**2/RESY(i)**2   *( YGOOD(I)*XGOOD(I) )
641    c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESX(i)**2
642    c$$$     +        *( XGOOD(I)*(1-YGOOD(I)) )
643    c$$$         print*,((XV(I)-XM(I))**2+(YV(I)-YM(I))**2)/RESY(i)**2
644    c$$$     +        *( (1-XGOOD(I))*YGOOD(I) )
645    c$$$         print*,XV(I),XM(I),XGOOD(I)
646    c$$$         print*,YV(I),YM(I),YGOOD(I)
647        ENDDO        ENDDO
648    c$$$      print*,'CHISQ ',chi2
649  *     ------------------------------------------------  *     ------------------------------------------------
650  *      *    
651  *     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 765  c$$$      ENDDO
765        RETURN        RETURN
766        END        END
767    
768    ******************************************************************************    
769    *    
770    *     routine to compute Likelihodd+Student and its derivatives
771    *    
772    *     (modified in respect to the previous one in order to include
773    *     single clusters. In this case the residual is evaluated by
774    *     calculating the distance between the track intersection and the
775    *     segment AB associated to the single cluster)
776    *
777    ******************************************************************************
778    
779          SUBROUTINE CHISQSTT(IFLAG,JFAIL)
780    
781          IMPLICIT DOUBLE PRECISION (A-H,O-Z)
782    
783          include 'commontracker.f' !tracker general common
784          include 'common_mini_2.f' !common for the tracking procedure
785    
786          LOGICAL TRKDEBUG,TRKVERBOSE
787          COMMON/TRKD/TRKDEBUG,TRKVERBOSE
788    
789          DIMENSION AL_P(5)
790          DIMENSION VECTEMP(5)
791    c$$$      DIMENSION U(5) ! BFGS
792    
793          DO I=1,5
794             AL_P(I)=AL(I)
795          ENDDO
796          JFAIL=0                   !error flag
797          CALL POSXYZ(AL_P,JFAIL)   !track intersection with tracking planes
798          IF(JFAIL.NE.0) THEN
799             IF(TRKVERBOSE)
800         $        PRINT *,'CHISQSTT ==> error from trk routine POSXYZ !!'
801             IFAIL=1
802             RETURN
803          ENDIF
804    
805          DO I=1,nplanes
806             DXDAL(I,1)=1.
807             DYDAL(I,1)=0.
808             DXDAL(I,2)=0.
809             DYDAL(I,2)=1.
810             COSTHE=DSQRT(1.-AL(3)**2)
811             IF(COSTHE.EQ.0.) THEN
812                IF(TRKVERBOSE)PRINT *,'=== WARNING ===> COSTHE=0'
813                IFAIL=1
814                RETURN
815             ENDIF
816             DXDAL(I,3)=(ZINI-ZM(I))*DCOS(AL(4))/COSTHE**3
817             DYDAL(I,3)=(ZINI-ZM(I))*DSIN(AL(4))/COSTHE**3
818             DXDAL(I,4)=-AL(3)*(ZINI-ZM(I))*DSIN(AL(4))/COSTHE
819             DYDAL(I,4)=AL(3)*(ZINI-ZM(I))*DCOS(AL(4))/COSTHE
820             IF(AL(5).NE.0.) THEN      
821                DXDAL(I,5)=
822         +           (XV(I)-(AL(1)+AL(3)/COSTHE*(ZINI-ZM(I))
823         +           *DCOS(AL(4))))/AL(5)
824                DYDAL(I,5)=
825         +           (YV(I)-(AL(2)+AL(3)/COSTHE*(ZINI-ZM(I))
826         +           *DSIN(AL(4))))/AL(5)
827             ELSE
828                DXDAL(I,5)=100.*( 0.25 *0.3*0.4*(0.01*(ZINI-ZM(I)))**2 )
829                DYDAL(I,5)=0.
830             ENDIF
831          ENDDO
832    
833          IF(IFLAG.EQ.0) THEN ! function calulation
834             CHI2=0.
835             DO I=1,nplanes
836                IF(XGOOD(I).EQ.1.AND.YGOOD(I).EQ.0)THEN !X-cl
837                   BETA = (XM_B(I)-XM_A(I))/(YM_B(I)-YM_A(I))
838                   ALFA = XM_A(I) - BETA * YM_A(I)
839                   YM(I) = ( YV(I) + BETA*XV(I) - BETA*ALFA )/(1+BETA**2)
840                   if(YM(I).lt.dmin1(YM_A(I),YM_B(I)))
841         $              YM(I)=dmin1(YM_A(I),YM_B(I))
842                   if(YM(I).gt.dmax1(YM_A(I),YM_B(I)))
843         $              YM(I)=dmax1(YM_A(I),YM_B(I))
844                   XM(I) = ALFA + BETA * YM(I) !<<<< measured coordinates
845                ELSEIF(XGOOD(I).EQ.0.AND.YGOOD(I).EQ.1)THEN !Y-cl
846                   BETA = (YM_B(I)-YM_A(I))/(XM_B(I)-XM_A(I))
847                   ALFA = YM_A(I) - BETA * XM_A(I)
848                   XM(I) = ( XV(I) + BETA*YV(I) - BETA*ALFA )/(1+BETA**2)
849                   if(XM(I).lt.dmin1(XM_A(I),XM_B(I)))
850         $              XM(I)=dmin1(XM_A(I),XM_B(I))
851                   if(XM(I).gt.dmax1(XM_A(I),XM_B(I)))
852         $              XM(I)=dmax1(XM_A(I),XM_B(I))
853                   YM(I) = ALFA + BETA * XM(I) !<<<< measured coordinates
854                ENDIF
855                TERMX = DLOG( (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)/
856         $           (TAILX(I)*RESX(I)**2) )
857                TERMY = DLOG( (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)/
858         $           (TAILY(I)*RESY(I)**2) )
859                CHI2=CHI2
860         $           +(TAILX(I)+1.0)*TERMX   *( XGOOD(I) )
861         $           +(TAILY(I)+1.0)*TERMY   *( YGOOD(I) )
862             ENDDO
863          ENDIF
864    
865          IF(IFLAG.EQ.1) THEN ! derivative calulation
866             DO I=1,5
867                CHI2DOLD(I)=CHI2D(I)
868             ENDDO
869             DO J=1,5
870                CHI2D(J)=0.      
871                DO I=1,nplanes
872                   CHI2D(J)=CHI2D(J)
873         $              +2.*(TAILX(I)+1.0)*(XV(I)-XM(I))/
874         $              (TAILX(I)*RESX(I)**2+(XV(I)-XM(I))**2)*
875         $              DXDAL(I,J)   *XGOOD(I)
876         $              +2.*(TAILY(I)+1.0)*(YV(I)-YM(I))/
877         $              (TAILY(I)*RESY(I)**2+(YV(I)-YM(I))**2)*
878         $              DYDAL(I,J)   *YGOOD(I)
879                ENDDO
880             ENDDO
881             DO K=1,5
882                VECTEMP(K)=0.
883                DO M=1,5
884                   VECTEMP(K) = VECTEMP(K) +
885         $              COV(K,M)/2.*(CHI2D(M)-CHI2DOLD(M))
886                ENDDO
887             ENDDO
888             DOWN1 = 0.
889             DO K=1,5
890                DOWN1 = DOWN1 + DAL(K)*(CHI2D(K)-CHI2DOLD(K))
891             ENDDO
892             IF(DOWN1.EQ.0.) THEN
893                PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN1 = 0'
894                IFAIL=1
895                RETURN
896             ENDIF
897             DOWN2 = 0.
898             DO K=1,5
899                DO M=1,5
900                   DOWN2 = DOWN2 + (CHI2D(K)-CHI2DOLD(K))*VECTEMP(K)
901                ENDDO
902             ENDDO
903             IF(DOWN2.EQ.0.) THEN
904                PRINT*,'WARNING IN MATRIX CALULATION (STUDENT), DOWN2 = 0'
905                IFAIL=1
906                RETURN
907             ENDIF
908    c$$$         DO K=1,5 ! BFGS
909    c$$$            U(K) = DAL(K)/DOWN1 - VECTEMP(K)/DOWN2
910    c$$$         ENDDO
911             DO I=1,5
912                DO J=1,5
913                   CHI2DD(I,J) = COV(I,J)/2.
914         $              +DAL(I)*DAL(J)/DOWN1
915         $              -VECTEMP(I)*VECTEMP(J)/DOWN2
916    c$$$     $              +DOWN2*U(I)*U(J) ! BFGS
917                ENDDO
918             ENDDO
919          ENDIF
920    
921          RETURN
922          END
923          
924  *****************************************************************  *****************************************************************
925  *  *
926  *     Routine to compute the track intersection points  *     Routine to compute the track intersection points
# Line 625  c      COMMON/TRKD/TRKVERBOSE Line 959  c      COMMON/TRKD/TRKVERBOSE
959  c        c      
960        DIMENSION AL_P(5)        DIMENSION AL_P(5)
961  *      *    
962        DO I=1,nplanes  cpp      DO I=1,nplanes
963           ZV(I)=ZM(I)            !  cpp         ZV(I)=ZM(I)            !
964        ENDDO        cpp      ENDDO      
965  *      *    
966  *     set parameters for GRKUTA  *     set parameters for GRKUTA
967  *      *    
# Line 641  c       Line 975  c      
975        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)        VOUT(6)=-1.*DSQRT(1.-AL_P(3)**2)
976        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))
977        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8        IF(AL_P(5).EQ.0.) VOUT(7)=1.E8
978    
979    c$$$      print*,'POSXY (prima) ',vout
980    
981        DO I=1,nplanes        DO I=1,nplanes
982           step=vout(3)-zv(i)  cpp         step=vout(3)-zv(i)
983             step=vout(3)-zm(i)
984   10      DO J=1,7   10      DO J=1,7
985              VECT(J)=VOUT(J)              VECT(J)=VOUT(J)
986              VECTINI(J)=VOUT(J)              VECTINI(J)=VOUT(J)
# Line 674  c$$$            if(.TRUE.)print*,'step', Line 1012  c$$$            if(.TRUE.)print*,'step',
1012              GOTO 11              GOTO 11
1013           ENDIF           ENDIF
1014    
1015    
1016  *     -----------------------------------------------  *     -----------------------------------------------
1017  *        evaluate track coordinates  *        evaluate track coordinates
1018   100     XV(I)=VOUT(1)   100     XV(I)=VOUT(1)
# Line 683  c$$$            if(.TRUE.)print*,'step', Line 1022  c$$$            if(.TRUE.)print*,'step',
1022           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)           AYV(I)=DATAN(VOUT(5)/VOUT(6))*180./ACOS(-1.)
1023  *     -----------------------------------------------  *     -----------------------------------------------
1024    
1025             IF(TRACKMODE.EQ.1) THEN
1026    *     -----------------------------------------------
1027    *        change of energy by bremsstrahlung for electrons
1028                VOUT(7) = VOUT(7) * 0.997 !0.9968
1029    *     -----------------------------------------------
1030             ENDIF
1031    
1032        ENDDO        ENDDO
1033    
1034    c$$$      print*,'POSXY (dopo) ',vout
1035    
1036    
1037        RETURN        RETURN
1038        END        END
1039    
# Line 725  c         ZM_B(IP) = 0 Line 1074  c         ZM_B(IP) = 0
1074           RESY(IP) = 1000.       !12.d-4           RESY(IP) = 1000.       !12.d-4
1075           XGOOD(IP) = 0           XGOOD(IP) = 0
1076           YGOOD(IP) = 0           YGOOD(IP) = 0
1077             DEDXTRK_X(IP) = 0
1078             DEDXTRK_Y(IP) = 0
1079             AXV(IP) = 0
1080             AYV(IP) = 0
1081             XV(IP) = -100
1082             YV(IP) = -100
1083        enddo        enddo
1084    
1085        return        return
# Line 795  c      IMPLICIT DOUBLE PRECISION (A-H,O- Line 1150  c      IMPLICIT DOUBLE PRECISION (A-H,O-
1150              ZP(NP)=ZM(I)              ZP(NP)=ZM(I)
1151           ENDIF           ENDIF
1152        ENDDO        ENDDO
1153          IFLAG=0                   !no debug mode
1154        CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)        CALL TRICIRCLE(NP,XP,ZP,AP,RP,CHI,XC,ZC,RADIUS,IFLAG)
1155  c      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI  
1156    c$$$      print*,' circle: ',XC,ZC,RADIUS,' --- ',CHI,IFLAG
1157    c$$$      print*,' XP ',(xp(i),i=1,np)
1158    c$$$      print*,' ZP ',(zp(i),i=1,np)
1159    c$$$      print*,' AP ',(ap(i),i=1,np)
1160    c$$$      print*,' XP ',(rp(i),i=1,np)
1161    
1162        IF(IFLAG.NE.0)GOTO 10 !straigth fit        IF(IFLAG.NE.0)GOTO 10 !straigth fit
1163    c      if(CHI.gt.100)GOTO 10 !straigth fit
1164        ARG = RADIUS**2-(ZINI-ZC)**2        ARG = RADIUS**2-(ZINI-ZC)**2
1165        IF(ARG.LT.0)GOTO 10       !straigth fit        IF(ARG.LT.0)GOTO 10       !straigth fit
1166        DC = SQRT(ARG)              DC = SQRT(ARG)      
# Line 806  c      print*,' circle: ',XC,ZC,RADIUS,' Line 1169  c      print*,' circle: ',XC,ZC,RADIUS,'
1169        AX = -(ZINI-ZC)/DC        AX = -(ZINI-ZC)/DC
1170        DEF=100./(RADIUS*0.3*0.43)        DEF=100./(RADIUS*0.3*0.43)
1171        IF(XC.GT.0)DEF=-DEF        IF(XC.GT.0)DEF=-DEF
1172          
1173    
1174    
1175          IF(ABS(X0).GT.30)THEN
1176    c$$$         PRINT*,'STRANGE GUESS: XC,ZC,R ',XC,ZC,RADIUS
1177    c$$$     $     ,' - CHI ',CHI,' - X0,AX,DEF ',X0,AX,DEF
1178             GOTO 10       !straigth fit
1179          ENDIF
1180        GOTO 20                   !guess is ok        GOTO 20                   !guess is ok
1181    
1182  *     ----------------------------------------  *     ----------------------------------------
# Line 847  c      print*,' circle: ',XC,ZC,RADIUS,' Line 1218  c      print*,' circle: ',XC,ZC,RADIUS,'
1218        AL(2) = Y0        AL(2) = Y0
1219        tath  = sqrt(AY**2+AX**2)        tath  = sqrt(AY**2+AX**2)
1220        AL(3) = tath/sqrt(1+tath**2)        AL(3) = tath/sqrt(1+tath**2)
1221        IF(AX.NE.0)THEN  c$$$      IF(AX.NE.0)THEN
1222           AL(4)= atan(AY/AX)  c$$$         AL(4)= atan(AY/AX)
1223        ELSE  c$$$      ELSE
1224           AL(4) = acos(-1.)/2  c$$$         AL(4) = acos(-1.)/2
1225           IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)  c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
1226    c$$$  ENDIF
1227    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
1228    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
1229    
1230    c$$$      AL(4) =  0.
1231    c$$$      IF(AX.NE.0.AND.AY.NE.0)THEN
1232    c$$$         AL(4)= atan(AY/AX)
1233    c$$$      ELSEIF(AY.EQ.0)THEN
1234    c$$$         AL(4) = 0.
1235    c$$$         IF(AX.LT.0)AL(4) = AL(4)+acos(-1.)        
1236    c$$$      ELSEIF(AX.EQ.0)THEN
1237    c$$$         AL(4) = acos(-1.)/2
1238    c$$$         IF(AY.LT.0)AL(4) = AL(4)+acos(-1.)
1239    c$$$      ENDIF
1240    c$$$      IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)
1241    c$$$      AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking ref.sys.
1242    
1243    c$$$      AL(4)=0.
1244    c$$$      IF( AX.NE.0.OR.AY.NE.0. ) THEN
1245    c$$$         AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
1246    c$$$         IF(AX.LT.0.) AL(4) = ACOS(-1.0)-AL(4)
1247    c$$$      ENDIF
1248    
1249          AL(4)=0.
1250          IF( AX.NE.0.OR.AY.NE.0. ) THEN
1251             AL(4) = ASIN(AY/SQRT(AX**2+AY**2))
1252             IF(AX.LT.0.AND.AY.GE.0) AL(4) = ACOS(-1.0)-AL(4)
1253             IF(AX.LT.0.AND.AY.LT.0) AL(4) = -ACOS(-1.0)-AL(4)
1254        ENDIF        ENDIF
1255        IF(AX.LT.0)AL(4)= acos(-1.)+ AL(4)        IF(AY.GT.0.) AL(4) = AL(4)-ACOS(-1.0)
1256        AL(4) = -acos(-1.) + AL(4) !from incidence direction to tracking rs        IF(AY.LE.0.) AL(4) = AL(4)+ACOS(-1.0)
1257    
1258        AL(5) = DEF        AL(5) = DEF
1259    
1260  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.17
  ViewVC Help
Powered by ViewVC 1.1.23